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Results from a recent quantum Monte Carlo (QMC) study (P.B. Chakraborty et al., Phys. Rev. 
B 70, 144411 (2004)) study of the LiHoF4 Ising magnetic material in an external transverse mag- 
netic field Bx show a discrepancy with the experimental results, even for small Bx where quantum 
fluctuations are small. This discrepancy persists asymptotically close to the classical ferromagnet to 
paramagnet phase transition. In this paper, we numerically reinvestigate the temperature T, versus 
transverse field phase diagram of LiHoF4 in the regime of weak Bx- In this regime, starting from 
an effective low-energy spin- 1/2 description of LiHoF4, we apply a cumulant expansion to derive 
an effective temperature-dependent classical Hamiltonian that incorporates perturbatively the small 
quantum fluctuations in the vicinity of the classical phase transition at Bx ~ 0. Via this effective 
classical Hamiltonian, we study the B^ — T phase diagram via classical Monte Carlo simulations. 
In particular, we investigate the influence on the phase diagram of various effects that may be at 
the source of the discrepancy between the previous QMC results and the experimental ones. For 
example, we consider two different ways of handling the long-range dipole-dipole interactions and 
explore how the Bx — T phase diagram is modified when using different microscopic crystal field 
Hamiltonians. The main conclusion of our work is that we fully reproduce the previous QMC re- 
sults at small Bx- Unfortunately, none of the modifications to the microscopic Hamiltonian that we 
explore are able to provide a Bx — T phase diagram compatible with the experiments in the small 
semi-classical Bx regime. 



I. INTRODUCTION 

A. Transverse Field Ising Model 

Phase transitions from order to disorder are most com- 
monly driven by thermal fluctuations. However, near 
absolute zero temperature, a system can, via quantum 
fluctuations associated with the Heisenberg uncertainty 
principle, undergo a quantum phase transition (QPT)ii^. 
The transverse field Ising model (TFIM) is perhaps the 
simplest model that exhibits a QPT— i^iii. This model was 
first proposed by de Gennes to describe proton tunneling 
in ferroelectric systems^. The Hamiltonian of the TFIM 
is given by 



H' 



TFIM 



(1) 



where erf (/Lt = x,y,z) are the Pauli matrices. Since 
af and af do not commute, a nonzero field F, trans- 
verse to the Ising z direction, causes quantum tunnel- 
ing between the spin-up and spin-down eigenstates of af , 
hence causing quantum spin fluctuations. These fluctu- 
ations decrease the critical temperature Tc at which the 
spins develop long-range order. In the simplest scenario, 
where Jy > 0, the ordered phase is ferromagnetii 
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At a critical field Fc, Tc vanishes, and a quantum phase 
transition between the quantum paramagnet (PM) and 
a long-range ordered ferromagnetic state occurs. The 
i/xFiM can be generalized by considering Jij as quenched 
(frozen) random interactions. Competing ferromagnetic 
Jij > and antiferromagnetic Jij < couplings gener- 
ates random frustration. For a three dimensional case. 



the system freezes into an (Ising) spin glass state at a spin 
glass critical temperature T^^. Similarly to the previ- 
ous example, Tg(T) decreases as F is increased until, at 
F = Fc, a quantum phase transition between a quantum 
paramagnet and a spin glass phase occurs. Extensive 
numerical studies have found the QPT between a quan- 
tum paramagnet and a spin glass phase^i^ii^ to be quite 
interesting due to the occurrence of Griffiths-McCoy sin- 
gularitiesiiii^. 



B. LiHo:,Yi :,F4 

The magnetic insulator LiHoF4, with a magnetic field 
Bx applied perpendicular to the Ising z direction of 
the Ho'^"'' magnetic moments, is a well known exam- 
ple of a physical realization of the transverse field Ising 
modeli^iii^ii^ii^. In LiHoF4 the predominant Jij interac- 
tion between the Ho'^^ ions is the long range interaction 
between magnetic dipoles which decays as l/'^fj, where 
rij is the distance between the i and j ions. The sign of 
Jij depends on the position of j respect to i- The exis- 
tence of a large crystal field anisotropy on the magnetic 
Ho'^"'' ionsi^ causes the system to behave as a classical 
Ising system with dipolar interactions for zero applied 
magnetic field - The reason is that the single ion crys- 
tal field ground state is an Ising doublet, meaning that 
the matrix elements of the raising and lowering angular 
momentum operator vanish within the space spanned 
by the two states of the doublet. The Ising direction 
is parallel to the c axis of the body centered tetrago- 
nal structure of LiHoF4. In zero applied magnetic field 



2 



B~c , the system is well described by a low-energy effective 
spin- 1/2 classical dipolar Ising modelilii^. Because the 
energy gap between the ground doublet and the first ex- 
cited singlet is fairly large compared to the Jy couplings, 
there is little quantum mechanical admixing between the 
ground doublet and the excited state induced by the in- 
teractionfiii. However, a nonzero admixes the ground 
doublet with the excited singlet and splits the ground 
doublet. It is this energy splitting which corresponds to 
the effective transverse field F in the TFIM description 
of LiHoF4 in nonzero iJ^J^ii^ii^. 

The Ho'^''' ions may be substituted (i.e. randomly di- 
luted) by non-magnetic yttrium (Y'^+) ions, with very 
little lattice distortion. This allows one to study the 
effects of disorder on LiHoxYi_j,F4 as an example of 
a diluted Ising model. Depending on the concentra- 
tion X of magnetic ions, the low temperature phase is 
either fcrromagneti a^^'^° or spin glass2ii22ii^. Interest- 
ingly, paradoxical behaviors are observed when a trans- 
verse magnetic field is applied to LiHo2;Yi_a;F4, with 
X < 1. In the ferromagnetic regime, (0.25 < a; < 1.0), 
when Ex = 0, a mean-field behavior Tc{x) oc x for 
the paramagnet to ferromagnet temperature transition 
is observed. However, in nonzero Bx, with increasing 
Bx, Tc{Bx) decreases faster than mean field theory pre- 
dicts^i. For Bx — 0, when LiHoxYi_i:F4 is diluted be- 
low X ~ 0.25, a conventional spin glass transition is ob- 
servedi^i^i^^. The signature of the spin glass transition 
is the divergence of the nonlinear magnetic susceptibil- 
ity X3 9't T^. However, surprisingly, XaiT) becomes 
less singular as Bx is increased from Bx = 0, suggest- 
ing that no quantum phase transition between a PM 
and a SG state exists as T — > 0^^^. Recently, theo- 
retical studiesiSii^i^Sii^ have suggested that for dipole- 
couplcd Ho'^+ in diluted LiHoj;Yi_2;F4, nonzero Bx gen- 
erates longitudinal (along the Ising z direction) random 
fields that couple to the magnetic moment and (i) lead to 
a faster decrease of Tc{Bx) in the ferromagnetic regime 
and (ii) destroy the paramagnet to spin glass transition 
in LiH0a;Yi_2,F4 samples that otherwise show a SG tran- 
sition when Bx = ^^'^^ . Recently, for the ferromagnetic 
regime, the influence of these induced random fields on 
the behavior of the linear magnetic susceptibility x in 
the presence of an external transverse magnetic field has 
been experimentally studied^S. When LiHoa;Yi_a;F4 is 
highly diluted (e.g. LiHoo.o45Yo.955F4), very interest- 
ing and peculiar behaviors are observed. AC suscepti- 
bility data show that the distribution of relaxation times 
narrows upon cooling below 300 mK ^^i^'^i'^^ . This be- 
havior is quite different from that observed in conven- 
tional spin glasses, where the distribution of relaxation 
times broadens upon approaching a spin glass transition 
at Tg > — This so-called antiglass behavior has been 
interpreted as evidence that the spin glass transition in 
LiHoj;Yi_3;F4 disappears at some nonzero Xc > O^i^. 
This is in contrast with theoretical arguments^^ which 
argue that, because of the long-ranged 1/r^ nature of 
dipolar interactions, classical dipolar Ising spin glasses 
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FIG. 1: The discrepancy between the experimental^^ phase 
diagram of LiHoF4 and quantum Monte Carlo (QMC) sim- 
ulations using stochastic series expansion for small Bx from 
Ref. [l^. The whole phase diagram is shown in the inset. 
At low temperature and high Bx, neglecting the large hyper- 
fine interaction A, generates a significant discrepancy between 
the experimental quantum critical point and the one obtained 
from simulation. However, at low Bx and close to the classical 
critical point, the hyperfine interaction is not a quantitatively 
important parameter. Other possibilities for the origin of this 
discrepancy have to be invoked in this regime. 

should have Tg{x) > for all a; > 0. However, recent nu- 
merical^i^i and experimental works^^ claim that a finite 
temperature paramagnetic to spin glass phase transition 
may not occur for x as large as Xc ~ 0.2. 

C. LiHoF4 as a TFIM 

In addition to the phenomena arising in the diluted 
regime of LiHo2,Yi_a;F4, the x = 1 regime also turns 
out to be interesting. There still exist problems for the 
pure LiHoF4, requiring the properties of this system in 
nonzero Bx to be rc-invcstigated more thoroughly. Per- 
haps surprisingly, it is just recently that the properties of 
LiHoF4 in a transverse external magnetic field have been 
studied in quantitative detail startin g fr om a truly micro- 
scopic spin Hamiltoniani^. In Ref. [l9|, which reported 
results from a quantum Monte Carlo (QMC) study us- 
ing the stochastic series expansion (SSE) technique^, a 
general qualitative agreement between the microscopic 
model and experimental datai^ was obtained. However, 
as illustrated in Fig. [1] there is significant quantitative 
discrepancy between the Monte Carlo results of Ref. [Toj 
and the experimental data of Ref [l^. In particular, 
the discrepancy between experiment and QMC results 
persists asymptotically close to the classical ferromag- 
netic to paramagnetic phase transition, where Bx/Tc and 
quantum fluctuations are perturbatively small. For very 
low temperatures and high Bx, it is crucial to consider 
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the hyperfine interaction in order to explain the behav- 
ior of the phase diagram close to the quantum critical 
poin^ii^i^. However, for very small Bx/Tc, the numer- 
ical results shown in Fig. [1] indicate that the effect of 
the hyperfine interaction is not important close to the 
classical transition at Tc- 

It was suggested in Ref. [l^ that this discrepancy 
between simulation and experiment, close to the clas- 
sical transition, may be related to some uncertainty in 
the crystal field parameters (CFP) used in the crystal 
field Hamiltonian, which enters in the TFIM description 
of LiHoF4, and which is simulated via QMC. Indeed a 
number of CFP sets obtained from different experimen- 
tal works, such as susceptibility measurements^^, neu- 
tron scatteringi^, and electron paramagnetic resonance 
experiments^, provide rather different values for the 
CFP. Specifically, different CFP would lead to different 
field {Bx) dependent effective coupling parameters in the 
TFIM description of LiIIoF4, which would result in dif- 
ferent Bx vs Tc phase diagrams. 

Yet, there are other factors of strictly computational 
nature which may be at the origin of the discrepancy 
illustrated in Fig. [TJ For example, because of the diffi- 
culties associated with dipolar interactions, calculations 
incorporating long-range dipolar interactions need to be 
performed quite carefully. Because of the long-range na- 
ture and angular dependence of dipolar interactions, the 
dipolar sum U{i) = —'^/N'}2,-{1 — Scos^ 9 ij)/rfj is con- 
ditionally convergen t'^^i'^^'^" , i.e the value of the sum de- 
pends on the shape of the external boundary of the sys- 
tem studied. Here, is the distance between site i and 
j, and 0ij is the angle between and the Ising spin axis. 

The conditional convergence of dipolar sums has been 
studied by Luttinger and Tiszai^. They performed the 
dipolar sum for a number of spin structures for systems 
with different external boundary shapes. For example, 
they considered an infinitely large system of dipoles on 
a body centered cubic lattice. They foimd that when 
the external boundary is spherical, the ground state is 
antiferromagnetic, while it is ferromagnetic for a needle- 
shaped sample. Later, Griffiths rigorously proved that 
for zero external field the free energy for a dipolar lat- 
tice system has to be independent of the sample shape in 
the thermodynamic limit^. The immediate consequence 
of Griffiths' theorem is that in zero external field, the 
net magnetization of the sample has to be zero. Other- 
wise, the field caused by the magnetic moments sitting 
on the boundary of the sample would couple to the dipo- 
lar moments of the sample, making the free energy shape 
dependent. Therefore, as a result of Griffiths' theorem^, 
domains must form in the sample, such the total mag- 
netization of the sample is zero in the thermodynamic 
limit. Griffiths' theorem is at variance with Luttinger 
and Tisza^S results because, in their work, the spin con- 
figurations were assumed uniform, and domain formation 
was neglected. This discussion emphasizes the complica- 
tion of studying systems with dipolar interactions and 
the caution which should be taken while dealing with 



such systems (e.g. the choice of the boundary geometry, 
boundary conditions and and the shape of the domains.) 
Finite size effects is another issue that needs to be han- 
dled quite carefully in systems where ions interact via 
long-range interactions. 

There are different ways to incorporate dipolar inter- 
actions in a computationally efhcient way. The method 
implemented in Ref. [l^ is the reaction field method^, 
which truncates the sum of the long-range interactions 
at the boundary of a sphere. The dipoles outside the 
sphere are treated in a mean-field fashion. Due to 
the semi mean-field nature of this method, the reac- 
tion field method overestimates the critical temperature. 
In the presence of quantum fluctuations, this overesti- 
mation is still at play and can possibly influence the 
Bx-Tc phase diagram as well. The Ewald summation 
mcthod^i^i^i^i^ is another method to treat the long- 
range dipolar interactions. In the Ewald summation 
method, a specified volume is periodically replicated. 
Then, by summing two convergent series effectively rep- 
resenting the dipolar interactions between magnetic mo- 
ments i and j, and all the periodically repeated im- 
ages of j, an effective dipolc-dipole interaction between 
two arbitrary magnetic moments i and j within the fi- 
nite size sample to be numerically simulated is derived. 
From a general perspective, it would appear quite worth- 
while to investigate the applicability and usefulness of the 
Ewald summation method to determine the low Bx vs Tc 
phase diagram of LiHoF4. Indeed, the Ewald summation 
method, unlike the reaction field one, is less prone to 
mean field over-estimations, and can be used as another 
methodology to probe the LiHoF4 problem via simula- 
tions^i. 

Another factor whose influence on the Bx — T phase 
diagram that should be studied is the nearest neighbor 
exchange interaction Jox in LiHoF4. The strength of 
•A;x, which is expected to be comparable to the dipo- 
lar interactions for a 4/ ion such as Ho'^^, is unknown. 
The strength can be determined such that the classical 
critical temperature matches the experimental value for 
Bx = 0. The estimated value of Jox is highly sensitive to 
the method used to handle the external boundaries and 
finite size effects in simulations, both of which have signif- 
icant effects when using the reaction field (RF) method, 
as already found in Ref. [l^ . 



D. Scope of the Paper 

The above discussion should make it clear that there 
are two rather distinct avenues to pursue in order to seek 
an explanation for the discrepancy between the experi- 
mentali^ Bx vs Tc phase diagram of LiHoF4 and the one 
obtained via QMGi^. One avenue, is that the current 
microscopic model is incomplete. As mentioned above, 
and suggested in Ref. [l^ , one possible source for this in- 
completeness may be an inaccurate set of CFP. Another 
possible source is that other interactions other than long- 
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range magnetic dipolar interactions and nearest neighbor 
exchange may be at plaj*^. Examples of other interac- 
tions include higher order multipole interactions and vir- 
tual phonon exchange^. The other avenue is related to 
the ensemble of computational pitfalls and insuing nu- 
merical errors that may arise when one deals with long 
range dipolar interactions through simulations. There- 
fore, before one delves into exploring a more complex 
microscopic Hamiltonian, there is a clear need to re- 
investigate the "simpler" problem that solely considers 
long-range dipole-dipole interactions and nearest neigh- 
bor exchange. 

In this work we aim to scrutinize the individual role 
of each of the computational issues as potential culprits 
for the discrepancy observed in Fig. [T] Because QMC 
and experiment do not match at B^jT^ — + 0, we have 
developed a tool that allow us to achieve the goal in an 
efficient and computationally simple way. Since this dis- 
crepancy appears at low enough near the classical Tc, 
where quantum fluctuations are perturbatively small, we 
can expand the partition function Z in terms of the trans- 
verse magnetic field B^ , and recast the partition function 
as a sum over strictly classical states, using a new ef- 
fective, albeit temperature dependent, classical Hamilto- 
nian Hcff{T). In Hcs(T), the quantum effects are incor- 
porated perturbatively, giving us the ability to calculate 
all thermodynamical quantities in presence of small quan- 
tum fluctuations within a classical Monte Carlo method. 
Therefore classical Monte Carlo simulations can be easily 
performed using Hag{T) in a very simple way, without 
the need to perform complicated QM C^^'^^ simulations 
when interested in a regime with weak quantum fluctua- 
tions^. Therefore, we can focus on the region close to the 
classical transition and investigate the different possible 
origins of the discrepancy in detail. 

In summary, (i) the complexity of the QMC SSE 
method, (ii) the problematic conditional convergence of 
dipolar lattice sums, (iii) the question of controlled finite 
size effects and its role on the consistent determination of 
the nearest-neighbor exchange Jox, and (iv) the possible 
sensitivity of the Tc{Bx) dependence on the choice of the 
CFP altogether warrant a new numerical investigation of 
the Tc{Bx) phase diagram in the LiHoF4 transverse field 
Ising material. Below, we will show that either fortu- 
nately or unfortunately, depending on one's disposition, 
the factors proposed in Section ITC] as the possible origins 
of the discrepancy between experiment and simulation 
(see Fig. [T]) are apparently not the issue. Therefore, the 
origin of the discrepancy remains unexplained. However, 
the pcrturbative cumulant Monte Carlo tool that we have 
devised can be used effectively to search for the cause of 
discrepancy. Without it, the discovery of the irrelevance 
of the above factors through a classical Monte Carlo sim- 
ulation would have been a more CPU time consuming 
burden. Ultimately, the same tool can also be used to 
explore the role of the small B^ when x ^ 9 20,27,28,29 ^ 
Indeed, constructing the whole x-Tc{B^) phase diagram 
in the "small i?^" vicinity of the classical x-Tc phase di- 



agram by performing solely classical Monte Carlo was 
an original key motivation for the development of the 
method presented in this paper. 

The rest of the paper is organized as follows. In Sec. 
II, we review the crystal structure and the physical prop- 
erties of LiHoF4 in a transverse field B^ and the effect 
of the choice of crystal field potential on the magnetic 
low energy states. In Sec. Ill, we introduce the full 
microscopic Hamiltonian of LiHoF4. We discuss how, 
for low energies, an effective spin- 1/2 Hamiltonian for 
LiHoF4 can be constructed, and explain how one can 
picture LiHoF4 in nonzero B^ as a dipolar TFIM. We 
then discuss how a semiclassical effective Hamiltonian 
is derived from the TFIM Hamiltonian by incorporating 
the transverse field term perturbatively via a cumulant 
expansion. In Sec. IV, we employ the semiclassical ef- 
fective Hamiltonian obtained in the previous section in 
classical Monte Carlo simulations for small B^- We dis- 
cuss the results obtained using either the reaction field 
or Ewald summation method for the long-range dipole 
interactions. We discuss how J^x is estimated and inves- 
tigate the sensitivity of the determined value upon the 
choice of the numerical method. Finally, we compare 
the Bx-Tc phase diagrams originating from two different 
sets of crystal field parameters. Section V summarizes 
our results. The paper also contains three appendices. 
Appendix A discusses details pertaining to the crystal 
field Hamiltonian. Appendix B gives some of the inter- 
mediate steps needed to construct the effective classical 
Hamiltonian T-lcs{T). Finally, Appendix C give the for- 
mulae needed to calculate physical thermodynamic quan- 
tities when doing classical Monte Carlo simulations with 

Hefi(r). 



II. STRUCTURE AND CRYSTAL FIELD 

The magnetic material LiHoF4 undergoes a second- 
order phase transition from a paramagnetic to a ferro- 
magnetic state at a critical temperature of 1.53 K— i^^. 
The critical temperature can be reduced by applying 
a magnetic field B^ transverse to the Ising easy-axis 
direction. The magnetic field induces quantum fluctua- 
tions such that beyond a critical field of B^ w 4.9 Tcsla, 
the system displays a quantum phase transition from a 
ferromagnetic state to a quantum paramagnetic state 
at zero temperatursi^. The magnetic properties of 
LiHoF4 are due to Ho'^+ rare earth magnetic ions. 
The electronic ground state of Ho^+ is 4/^°, which 
gives small exchange couplin g^^i^°'^^ , such that the 
predominant magnetic interaction between the Ho'^"'" 
ions are long-range magnetic dipole-dipole interactions. 
Hund's rules dictate that the total angular momentum 
of a free ion }io^+, J = 8 (L = 6 and S = 2) and the 
electronic ground state configuration is ^1$- LiHoF4 is 
a compound with space-group (/4i/a) and lattice 
parameters a = b = 5.175A, c ~ 10.75A, and has 4 Ho'^+ 
ions per unit cell positioned at (0,0,1/2), (0,1/2,3/4), 
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(1/2,1/2,0) and (1/2,0,1/4)^. The crystal has S4 
symmetry, which means the lattice is invariant with 
respect to a ^ rotation about the z axis and reflection 
with respect to the x — y plane. 

In the crystal structure, the Ho'^+ ions are surrounded 
by F~ ions, which create a strong crystal electric field 
with Si symmetry. This crystal field lifts the 17- fold de- 
generacy of the ^/g configuration giving a non-Kramers 
ground state doublet. The next excited state is a singlet 
with an energy gap of « 11 K above the ground state 
doubletji^'^i^'^. The crystal field Hamiltonian and the 
crystal field parametrization is discussed in more detail 
in Appendix [X] Holmium is an isotopically pure element 
with nuclear spin I = 7/2, which is coupled to the elec- 
tronic spin J via the hyperfine contact interaction Al ■ J, 
where A w 39 mK^^^. 




a = b = 5.116 k 
c = 10.75 A 



FIG. 2: The crystal structure of LiHoF4. NN identifies the 
first nearest neighbors and NNN identifies the next nearest 
neighbors 



III. EFFECTIVE THEORY OF LiHoF4 FOR 
THE LOW Bx/Tc REGIME 

In this section we derive an effective model suitable for 
describing LiIIoF4 in a small transverse magnetic field 
regime, where Bx/T^ —f {Tc is the critical temperature 
when Bx = 0). The simplicity gained using an effec- 
tive theory gives us the ability to capture the essential 
physics, and to easily reinvestigate the influence of the 
different parameters affecting the behavior of the phase 
diagram of LiIIoF4 in the B^/Tc regime. We de- 



rive the required effective model in two steps. Firstly, 
in LiIIoF4, in the temperature range that we are inter- 
ested in, which is close or below Tc{Bx = 0) = 1.53 K, the 
high energy scales are well separated from the low energy 
sector. The energy scale for dipolar interactions between 
nearest-neighbor Ho'^+ ions is about 0.31 K. This is much 
smaller than the energy gap between the two first lowest 
single ion energy states and the next higher crystal field 
states (> 11 K). In this case, one can neglect the higher 
energy states and reduce the full Hamiltonian Hilbert 
space to a smaller subspace spanned by the two lowest 
energy states. This enables us to deduce a low energy ef- 
fective spin-i Hamiltonian for LiHoF4. Secondly, we de- 
rive a semi-classical effective Hamiltonian from this low 
energy spin-^ Hamiltonian by incorporating the trans- 
verse field term pcrturbativcly via a cumulant expansion. 
We can then perform a simple classical Monte Carlo us- 
ing this semi-classical effective Hamiltonian to investigate 
the small Bx/Tc regime. 

A. Effective Spin-i Hamiltonian 

As mentioned in the previous section, there are three 
type of interactions that play a role in the magnetic prop- 
erties of LiHoF4. The main interaction is the long-range 
dipole-dipole interaction between the Ho'^^ magnetic ions 
denoted by 

H^^^ = \{g^^,^fY.T.L^P'^h (2) 

where iJ,i'=x,y,z and is the total angular mo- 
mentum of Ho'^^ ion i. L'^^ is the magnetic 
dipole interaction written in the form L^^ = 
[<5^''|ryf-3(ry)^(r,y)1 /|ryf , where ry is the dis- 
tance between ion i and j. gh = 1-25 is the Lande g-factor 
of free Ho'^^ and fiB = 0.6717 K/T is the Bohr magne- 
ton. The dipolar interaction is complemented by a short 
range nearest-neighbor Heiscnberg exchange interaction 

i,NN 

where NN denotes the nearest neighbors of site i. 
This exchange interaction is considered to be weak and 
isotropioi^i^. The third interaction is the hyperfine cou- 
pling between the electronic and nuclear magnetic mo- 
ments 

i/hyp=A^(I,.J,) . (4) 

i 

The hyperfine constant ^ « 39 mK is anomalously 
large in Ho'^^-bascd material o^'^'^^i'^^ . Thus, the complete 
Hamiltonian is written as 

i i 
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(Tesla) 



FIG. 3: The energy splitting of the ground state doublet, 
A{Bx) = Ei3{Bx) — Ea{Bx), in LiHoF4 as a function of the 
transverse magne tic field. The crystal field Vc was obtained 
from Refs. ^\^. For more details on the crystal field and 
crystal field parametrization (see Appendix IX)) . 



Projecting the single ion Hamiltonian of Eq. ([6]) in this 
two-dimensional subspace for an arbitrary ion z, we get 



(8) 



where E{B,) ^ \{E^{Bx) + EpiB^)) and A(B,) = 
Ef3{Bx) ~~ Ea[Bx)- The energy difference between the 
two lowest states caused by the transverse magnetic field 
Bx can already be interpreted as an effective transverse 
field F = A{Bx)/2 acting on ScS=l de grees of freedom 
at each site. The dependence of A(i?^) on the magnetic 
transverse field B^ is plotted in Fig. [3l 

Since we are henceforth working in a two-dimensional 
subspace for each ion i, we can write the interactions 
between 3f and J!f in terms of effective interactions be- 
tween Pauli matrices. Indeed, any operator acting in a 
two-dimensional space can be written as a linear combi- 
nation of (T^ Pauli matrices plus the unit matrix cr*^ = 1. 
In order to express Jf in terms of erf, we project Jf in 



the subspace spanned by 
write the operator as 



t) and II). Specifically, we 



The first two terms are single ion interactions, where Vc 
describes the strong crystal field interactions discussed in 
Section Inland Appendix lAl The second term is the Zee- 
man interaction. Henceforth, we ignore i^hyp since our 
goal, as explained in the Introduction, is to investigate 
the small B^ and smaU (Tc(0) - Tc{Bx)) /Tc(0) regime 
where, as already suggested by the results of Ref. and 
shown in Fig.[Tl the hyperfine interaction effects are neg- 
ligible. The first two single-site (non-interacting) terms 
in H , denoted as 



single- 



-site = Vb(J) - .gLMB^ajJ" 



(6) 



can be easily numerically diagonalized for arbitrary 
transverse field Bj^. \a{Bx)) and \l3{Bx)) are the two 
lowest states of the single ion Hamiltonian © for a given 
Bx- Their corresponding energies are denoted by Ea{Bx) 
and Ei3{Bx). At Bx = these two states form a doublet, 
but Bx ^ lifts the degeneracy. The Ising subspace 1 1) 
and 1 1) are chosen by performing a unitary rotation on 
the \a{Bx)) and \l3{Bx)) states : 



T) 



1 

1 

71 



(|q.) -cxp(z0)|/3)) 



(7) 



The phase 9 is chosen such that the matrix elements of 
the operator between | f) and | j) is real and diagonal, 
giving for Jf, Jf = Czzcrf ■ Since the first excited state, 
\^{Bx)), above \a{Bx)) and \(3{Bx)), is at an energy at 
least seven times higher than kBTc{Bx), and is repelled 
for all Bx from the \a{Bx)) and \P{Bx)) set (see Fig. 1 of 
Ref. [l^l), we henceforth neglect all excited crystal field 
states and work in a reduced Hilbert space spanned solely 
by \a{Bx)) and \l3{Bx)), or equivalently by 1 1) and | i). 



where 





I'—x^y 


z 




= ^[(TIJIT)- 


(ilJli)] , 




= ^[(T|J1T) + 




C ^x 


= ^[(T|Jni) + 


(ilJ^T)] and 




= ^KTlJ^i)- 


aiJiT)] . 



(9) 



Based on the crystal field parameters of Refs. [313, 
the evolution of the various parameters C^^ and C^jo as 
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FIG. 4: The evolution of the C^v parameters using the crystal 
field Vc from Refs. [l5lll9| . In the inset one can see that Cxy ~ 
Cyo. Coefficients that are not plotted are zero. 
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a function of is plotted in Fig. [4] We see that Czz is 
the largest term compared to all the other C^^'s. 

For the Hamiltonian in Eq. the operators are 
substituted by their two dimensional representations in- 
troduced in Eq. ([9]). This leads to a complicated Hamil- 
tonian that acts within the Ising subspace of | "f) and 
I I). The projection generates various kinds of interac- 



tions among the effective Scff = 5 spins. Via Eq. ([7]), a 
specific rotated subspace was chosen, such that C^^ = 
{fi — x,y,0; cr" = 1). As shown in the inset of Fig. HI 
Cxy, Cyx, and Cyo are very small, so the interacting terms 
containing these coefficients can be neglected. Therefore, 
neglecting these terms, we obtain 



|JoxE^'m(^-) E ^f'^NN + {9L^^BrCzz{B,)CMB.)Y.^^at 

fj, i.NN i^j 



E 



(10) 



When the external magnetic field B^ is zero, only 
Czz(O) 7^ and all the other C^i, and C^o vanish. Hence, 
in absence of an external magnetic field, the system can 
be described by a simple classical dipolar Ising modeU^. 
Fortunately, a number of interaction terms are zero or 
can be neglected with respect to the leading Ising inter- 
action, which is proportional to C^ziBx)J2i^j ^ij'^i'^j- 
As we can see from Eq. ([TU]). for pure LiHoF4, an effec- 
tive crfcrj and u^a-- pair- wise interactions as well as a 
linear transverse field along the x direction are induced 
in the presence of an external magnetic field. As sug- 
gested by Fig. [5l and already assumed in Ref. [l^ , we 
expect the quantum fluctuations induced by these terms 
via either dipolar or exchange coupling, to be quite small 
and negligible compared to the quantum fluctuations in- 
duced by A{Bx). For the pure (disorder free) LiHoF4, the 
invariance of the dipolar interactions under lattice mir- 
ror symmetries forces J^j ^Ij =0. So the linear term 
with Czz{Bx)Cxo{Bx)Y,i^j Llf<7i vanishes. Consider- 
ing the Czz{Bx)C,^.x{Bx) j2i^j Llj^i^'j term, because of 
lattice mirror symmetry, one has 'Ylii^j ■^tj'^i (^j) ~ 
0, therefore this term can only contribute via high 
order fluctuation effects beyond the vanishing mcan- 
field contribution. Since ^'""^^^j < 1, we expect the 
(second order) fluctuation contribution effects from the 
above cfcj term to be small. Hence we neglect the 
Cz.iBx)CxxiBx)J:^^JL|^<y!a^ term in the SeS = i ef- 
fectivc Hamiltonian i?spin-i/2- We should emphasize 
that for diluted LiHo3:Yi_2^F4, since the lattice mirror 
symmetries are broken, the two latter terms, proportional 
to 'Ei^j Llj^t and Y.^^j ^f/^rf (^^J ): can no longer be 
neglected^^. Indeed, these are the terms responsible 




0' ' ' ' ' ' 
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(Tesla) 

FIG. 5: The ratio of the typical value of terms neglected 
in Hamiltonia n (I lip respect to A, using the crystal field Vc 
from Refs. [islligj and the dipolar sum is performed for a long 
cylindrical sample. 



for the generation of the longitudinal random fields in 
LiHo2;Yi_a;F4 when subject to nonzero Bj^^S^J^^ as dis- 
cussed in the Introduction. 

Hence, the spin-i Hamiltonian in Eq. (jTU]) can be 
further simplified to a familiar looking transverse field 
Ising Hamiltonian with a dipolar and nearest-neighbor 
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exchange Ising interaction. 



1 



spin— 1 /2 



z z 
"NN 



:.NN 



To simplify the calculations, and in order to be consis- 
tent with the notation of Ref. [l^ as well as for further 
comparison between our simulation results and those of 
Ref. (loj , we lump the whole dependence in the trans- 
verse field term into a renormalization factor e{Bx) is 
defined as 



Czz{Bx) 



c^zio) ■ 

We renormalize the Hamiltonian as 

^spin-l/2 = [<^{Bx)]'^H 

with, according to Eq. pT|) . H is 



(12) 



(13) 



i,NN 



-.9lMbC,,(0)S,E< 



(14) 



where the renormalized effective transverse magnetic 
field Bx, is related to the real applied B^ via 



Bx 



AiBx) 



2,9LAiBC..(0) X [eiBx)Y 



(15) 



consistent with Ref. [l9[ . In discussing Monte Carlo sim- 
ulations below, we also define a renormalized tempera- 
ture, r, in conjunction with 7i, with T defined as 



T=[e{Bx)fT, 



(16) 



where T is the real physical temperature. 

All results presented in the Monte Carlo simulations 
section below were obtained by considering the renor- 
malized Hamiltonian (|14[) . and performing the simula- 
tions with respect to the renormalized T and Bx- Before 
presenting our Monte Carlo simulations of Eq. (fH|) as 
pertain to LiHoF4, we first discuss the technique we em- 
ployed to handle quantum fluctuations pcrturbatively for 
small Bx/Tc. 



B. Effective classical temperature-dependent 
Hamiltonian perturbation expansion 

In this section, with a focus on the simplified spin i 
Hamiltonian of Eq. (|14p . we aim to implement a cumu- 
lant pcrturbative Monte Carlo method for a spin i trans- 
verse Ising model^i^. For small quantum fluctuations. 



close to the classical critical temperature, we are able 
to derive an effective classical Hamiltonian analytically, 
where quantum fluctuations are incorporated perturba- 
tively. Using such effective perturbative Hamiltonian, we 
can then perform classical MC simulations. To set the 
stage, we first consider a general transverse field Ising 
Hamiltonian such as 



n 



1 



id 



j.NN 



(17) 



r is the transverse field in the x direction and ho de- 
notes an external longitudinal field along the z direction. 
For compactness, note that we passed from dipolar in- 
teractions denoted C'^^{0){ghfJ-B)^Lfj to £fj and from 
exchange interaction C^^(0)Jox to jTcx ( see Eq. (fT4|) ). 
The partition function Z for a system with Hamiltonian 
(HIl) is 



Z = Trace(e-'3^) 



(18) 



where Z is obtained by tracing over ^/j^'s which are, for 
example, direct product of erf eigenvectors (| t) and | J,)) 
and f3 = 1/kBT. We can write the Hamiltonian (|17|1 as 
H = Ho + Hi. Hq is the classical part of the Hamilto- 
nian, for which the "^i's are eigenvectors. Hi = — F J2i "'f 
is the quantum term, which does not commute with Ho- 
The existence of these two non-commuting terms in Ti. 
prevents us from applying classical Monte Carlo tech- 
niques directly to the system. We can derive an effective 
classical Hamiltonian as a functional of tjji, such that 



(19) 



Referring to the definition above in Eq. p9p . since the 
right hand side of Eq. (|19p is the matrix element with 
respect to {i/ji), H^s [ipi] is a functional depending only 
on the set of trf eigenvalues. The partition function can 
then be written as a classical partition function 



Z : 



(20) 



By finding an explicit expression for Hcs [ipi\j one can 
calculate the thermodynamical properties of the system 
described by Ti by performing classical Monte Carlo 
simulations using HcS instead of Ti.. 

To proceed, we write the matrix element {ip\e~^'^\ip) 
in terms of a cumulant expansion^ 



exp 



n>l 



(21) 
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To make the notation more compact, by we mean a 
typical eigenvector. Using Eq. (PT|) we can derive 
the effective Hamiltonian H^s [tl^i] perturbatively. The 
details of the derivation of ["(Aj] are presented in Ap- 
pendix |B] iJoff [V-'iJj is to order O(r^), given by 

i 

~Fo[2f3{h, + h„)]}. (22) 

In Eq. (im, hi is the total local field affecting the spin at 
site i caused by all the other spins, and which is 



Jc: 



NN 



(23) 



and ho is the external longitudinal field in the z direction. 
The functions Fo{x) and Fi{x) are defined as 



Foix) 
F,{x) 



_ cosh(a;) — 1 
_ sinh(a;) — x 



(24) 



In this effective Hamiltonian, the effect of quantum 
fluctuations is taken into account perturbatively to order 
0{pr'^ /[Hq]), where [Hq] denotes the order of magnitude 
of Ho, the classical part (first two terms) of Eq. (fT?]) . 
To obtain the thermodynamical properties of the sys- 
tem for small transverse fields wc can therefore perform 
a classical Monte-Carlo on Hctf as a classical counter- 
part of the real quantum mechanical Hamiltonian. Since 
we are interested in thermal averages we can calculate 
thermodynamical quantities by differentiating the parti- 
tion function, which is written in terms of H^s [ipi], with 
respect to ho, T or (3. The effective Hamiltonian has 
an explicit ho and /? dependence. For each true ther- 
modynamical quantum-mechanical quantity, we obtain a 
pseudo-operator counterpart. For example the pseudo- 
operators corresponding to (E), {AL), (M^), (M^), and 
(M^) are calculated in Appendix [Cl where E, and 
are the energy and magnetization operators along 
the z and x direction. (...) stands for the Boltzmann 
thermal average. 

Because of its perturbative nature in (/3r), this method 
is not reliable for large transverse fields or low tem- 
peratures. To illustrate the range of vahdity of this 
method we consider a simple one-dimensional nearest- 
neighbor transverse-field Ising-model Hamiltonian H = 
— JJ2i'^i'^i+i ~ with periodic boundary condi- 

tions. For a one-dimensional chain of 10 ions, wc arc 
able to calculate the exact total energy of the chain by 
exact diagonalization. To check our perturbative MC 
technique, we calculated the energy of the Ising chain 
as a function of temperature for a given transverse field. 
To make a comparison, we also performed a quantum 
Monte-Carlo (QMC) simulation on the system. In this 
QMC simulation, we used the Trotter-Suzuki^ formal- 
ism and applied a continuous time cluster algorithm sim- 
ilar to the one in Rcf . [5§| . In Fig. [HI for a quite large 



transverse field T/J = 1, we plot the average thermal 
energy as a function of temperature obtained from exact 
diagonalization, time cluster QM and "perturbative MC" 
using the effective perturbative Hamiltonian described 
above. This tests confirms the quantitative correctness 
of the perturbative Monte Carlo scheme at small /JF^/J. 
We also computed other thermodynamic quantities (e.g. 
(Mz), (Mx) ) and these also compared well with QMC 
and exact diagonalization results. 



-0.5 



E/J 
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Perturbative MC 
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FIG. 6; Energy as a function of temperature for a simple one 
dimensional nearest-neighbor Ising chain with a transverse 
field of r = J and = 10 spins and periodic boundary 
conditions. The energy is obtained by exact diagonalization 
of the Hamiltonian, a time-cluster QMC algorithm, and a 
classical Monte-Carlo algorithm of the perturbative effective 
Hamiltonian. 

Before we present our Monte Carlo results for LiHoF4, 
let us summarize what we have done so far. 

1. Since the spin-spin interactions and Tc{Bx) are 
small compared to the gap between the low-lying 
states \a{Bx)) and \(3{Bx)) with respect to the ex- 
cited crystal field state \^{Bx)), we can recast the 
full microscopic model of LiHoF4 in terms of an 
effective transverse field Ising model with effective 
spin-spin interactions and effective transverse field 
T(Bx) that depend on the real physical applied 
magnetic field B^. 

2. Since we are interested in a regime where B^/Tc is 
small, we can develop a perturbation expansion of 
the partition function in powers of B^/T and recast 
the thermal averages of real physical observables in 
terms of quantities that can be determined via a 
classical Monte Carlo simulation of a further effec- 
tive temperature-dependent classical Hamiltonian. 

Having shown that the perturbative cumulant MC can 
quantitatively describe the TFIM for small PT'^/[Ho], we 
proceed in the next section to describe how we use this 
method to study LiHoF4 at small transverse field B^, 
BjTc « 1. 
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IV. PERTURBATIVE MONTE CARLO STUDY 
OF LiHoF4 



and the external magnetic field acting on the domain is 
Bf^^, then the susceptibility x of the domain is 



In this section we report results from the perturbative 
Monte Carlo (MC) simulation to study the low trans- 
verse field Bx properties of LiHoF4, using the low field 
perturbative effective Hamiltonian in Eq. (|22|) and using 
Eq. for the definition of the local hi fields. As dis- 
cussed in the Introduction, our primary goal here is to 
check the quantum Monte Carlo results from stochastic 
series expansion of Ref . , and investigate the contrast- 
ing results with the transverse field B^ phase diagram of 
Ref. for small B^ (See Fig.[T]). Hence, we are indeed 
interested in LiIIoF4 in the case of asymptotically small 
Bx/Tc- The temperature we use in our simulations is the 
renormalized temperature defined in Eq. p6p . Regard- 
ing Eq. p4)) . the transverse field F used in the pertur- 
bative efl^ective Hamiltonian ((22|) is F = gLfJ'BCzz{0)l3x, 
where Bx is defined in Eq. p^ . For the local field hi, de- 
fined in Eq. (Eg), we have Lf^ = C^,,{0){gLm? Lf^' and 

In the following subsections, we first discuss the reac- 
tion field (RF) and the Ewald summation (ES) methods 
that we use to deal with the long range dipolar inter- 
actions, and discuss how the Monte Carlo results in the 
classical regime, where Bx ~ 0, are affected by the choice 
of the method we use. Next, we discuss the sensitivity of 
the Jox estimates at zero Bx to finite-size effects, bound- 
ary conditions and choice of the method to handle the 
dipolar lattice sum. We also consider the effect of differ- 
ent Jox on the phase digram, when Bx ^ and Bx/T is 
small. Finally, we investigate to what extent the final re- 
sults depend on the set of crystal field parameters chosen 
to describe the Ho'^"'" single ion properties. 



A. Reaction Field Method vs Ewald Summation 
Method 

Griffiths' theorem^ states that in the absence of an ex- 
ternal field the free energy for a dipolar lattice system has 
to be independent of the sample shape in the thermody- 
namical limit. Therefore, as an immediate consequence, 
in the absence of an external field, the net magnetization 
A4 of the sample has to be zero. Otherwise, for a uniform 
M. ^ 0, a, shape dependent demagnetization field would 
couple to the dipolar moments of the sample, making the 
free energy shape dependent. Here, the demagnetization 
field is the field originating from the magnetic moments 
sitting on the boundary of the sample. Hence, in the 
thermodynamic limit, domains form in order for the sys- 
tem to have a zero magnetization, M = 0. 

Experiments on LiHoF4 show that the results are 
shape independent, confirming Griffiths theorem and do- 
main formation^i^. There is evidence that in LiHoF4 
long needle-shaped domains form along the c axis^SiSi. 
If we assume that there is a uniform macroscopic bulk 
magnetization A4z within a long needle-shaped domain 



X = Mz/Bl 



(25) 



It should be noted that the macroscopic bulk magne- 
tization A^z, is given by A4z = 'T'OSlMb {J^') , where 
no = is the number of dipoles per unit of vol- 

ume and where a^c is the volume of the unit cell. Using 
= CzzCT^, the bulk magnetization Mz is related to the 
total moment of the effective Ising spins, Mz ~ 
in the S'cff=l/2 picture by 



iV a^c 



(26) 



where N is the total number of dipoles. 

Let us consider consider an imaginary macroscopic 
spherical cavity deep inside a needle-shaped domain. The 
magnetization inside the sphere should be equal to the 
uniform bulk magnetization of the long needle-shaped 
domain. Apart from the external magnetic field B^^*", 
spins enclosed in the sphere experience an additional field 
that originates from the spins on the outer boundary 
of the imaginary sphere embedded in the long needle- 
shaped domain. The magnetic surface charge density 
on the surface of the needle-shaped domain with uni- 
form magnetization A4z produces an internal magnetic 
field -Bneodic = 47rA^z- Meanwhile, the magnetic surface 
charge density on the surface of the uniformly magnetized 
sphere with magnetization of A^z induces a (demagneti- 
zation) magnetic field ^Aiz inside the sphere that is in 
the opposite direction to the applied field and to iJneedie- 
Therefore, the total field B^^^ inside the spherical cavity 



(27) 



Aiz is uniform for a bulk sample. Now, instead of con- 
sidering a whole needle-shaped bulk, we can also study 
an isolated spherical sample which an effective B^^^ field 
is applied to it. If we substitute Bf^^ with Mz/x ^^d 
^sph ^j^j^ Mz/Xsph, where Xsph is the susceptibility of 
the spherical domain, then we can write x as a function 

of Xsph 



Xsph 



3 Xsph 



(28) 



If Xsph is obtained via some calculation procedure for 
a spherical sample, one can use Eq. (|28p to determine 
the macroscopic susceptibility of the bulk sample within 
which the sphere is embedded. Specifically, simulations 
can be performed on a finite size sphere, and the effect of 
the macroscopic bulk surrounding the sphere is incorpo- 
rated in a mean-field manner by considering an effective 
field _B|P'^ interacting with the spins inside the spheri- 
cal sample. Using this method, called the reaction field 
(RF) method, Chakraborty et al. calculated the finite 
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size sphere susceptibility Xsph by using the stochastic se- 
ries expansion quantum Monte-Carlo metho d^^i'^^ . They 
considered an A'^ spin system enclosed by a sphere, where 
the susceptibility of the sphere is obtained from the spin- 
spin correlation. Referring to Eq. (|28p . the paramagnetic 
to ferromagnetic transition (criticality) within the macro- 
scopic long needle-shaped domain occurs at the temper- 
ature for which Xsph = ^ occurs for a spherical sample. 
It should be noted that this criteria is derived for macro- 
scopic systems in the thermodynamic limit. Therefore, 
as discussed in Ref. [i^, because of the fluctuation of 
magnetic moments on the boundary of a finite size sur- 
face, quantities such as specific heat and susceptibility 
obtained via the RF method, are quite sensitive to finite 
size effects. 

The Ewald summation (ES) method^i^i^ is an al- 
ternative approach used to obtain reliable quantitative 
results for describing real dipolar materials in a periodic 
boundary condition (PBC)^i^. In the ES method, in 
order to treat long-range dipolar interactions with PBC, 
the system is modeled by replicating the simulation cell 
of linear size L into a large array of image copies. The 
ES method generates an effective dipolc-dipole interac- 
tion ^ Lg^(ry)/i,f /ij between two arbitrary magnetic 
moments, /x^ and (ij within the simulation cell. Here, 
Mi = ShfJ-B^i, fi, i'=x,y, z, and Tj — where is the 
position of moment i. This is done by periodically repli- 
cating the simulation cell with a volume of Hq = L^a^c 
and summing convcrgcntly the interactions between the 
real spins i and j in the specified volume of the simulation 
cell of size L, and all the periodically repeated images of 
j as 



(29) 



where n = {nxLa,nyLa,nzLc) with rix, ny, riz integers. 



3(ri,)'^(r,,r] /|r„f are 



dipolar couplings, which can be written in a more com- 



pact form as L'""(ry) = Vf V^'lr,, h^. Therefore 



i:;^(r.,) = VrV^E|r,,+n|-i 



(30) 



The sum jr^ + n|^^ is calculated using the Ewald 
method, such that the sum contain a real space sum plus 
a reciprocal space sum minus a self tcr m^^i^^i^^ 



El 



Eerfc(K|ry + n| 
Ir.-. -I- nl 



e ' cos (k • Tij) 
(31) 



Here erfc(x) = {2/^/tt) x exp —t^dt and k denotes the 
reciprocal vectors of the simulation cell. The convergence 
factor K is chosen such that the real space sum and the re- 
ciprocal space sum converge about equally rapidly^i^i^. 



The simulation cell and all its replicated images are em- 
bedded altogether in a continuous medium. Additionally, 
each spin experiences a demagnetization field, which is 
originating from the magnetic moments on the boundary 
of the system^. This boundary contribution depends 
on the shape of the boundary of the macroscopic sam- 
ple that we are interested in modeling, i.e. for a long 
needle-shaped sample the demagnetization field correc- 
tion to the ES representation of the dipole-dipole inter- 
actions is zero — . However, for a bulk spherical sam- 
ple, the magnetic polarization of the magnetic moments 
on the boundary of the sphere induces a demagnetiza- 
tion field proportional to the magnetization of the sam- 
ple Ai = j^^i fJ-i: which creates an additional effective 
field acting on the the magnetic moments. The net effect 
results in an extra effective interaction 



2fi' + 1 Qq 



(32) 



between magnetic moments /x^ and fij to be incorporated 
in the simulation^. In practive, the term in Eq. (j32p is 
merely added to ^eff ('"u ) ^'i- which itself is cal- 
culated via the ES expression of Eq. pT|) . Here, L is the 
linear system size, fii = gLA*BJi, and /i' is the magnetic 
permeability of the surrounding continuum. For a sam- 
ple surrounded by vacuum /i' = 1—. This interactions is 
added to the effective dipolar interaction between spins 
i and j, derived by the ES technique^. 

As a result, within the ES method, each spin interacts 
with all the "real" spins in the specified simulation cell 
of linear size L, and with all its replicated periodic im- 
ages. Therefore, one would expect the system to behave 
more like a macroscopic system than in the RF method. 
However, there are still some finite size effects due to 
the artifact of having a periodic sequence of cells of fi- 
nite size L. Once an effective dipole-dipole interaction 
between spins i and j within the simulation cell has been 
derived via the ES technique, one can perform Monte 
Carlo simulations using the standard Metropolis algo- 
rithm. Xu et al^ used this ES technique to simulate 
long-range dipolar Ising interactions for both the body- 
centered cubic (BCC) and body centered tetragonal lat- 
tices in zero applied field. In a more recent workSl^ the 
ES technique was implemented in a Monte Carlo simula- 
tion study of LiHoa;Yi_a;F4 in zero applied field. In the 
next subsection we discuss the results of MC simulations 
using the cumulant perturbative method. In our simula- 
tions, we incorporate the long-range dipolar interactions 
using both the RF method as discussed in Ref. [l^ and 
the ES method. The influence of each method on the 
MC results is investigated in some detail. 



B. Perturbative Monte Carlo Simulations Results 

In this subsection we describe the Monte-Carlo re- 
sults obtained using the effective perturbative Hamilto- 
nian (j22p and which employ different ways to handle the 
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dipolar lattice sums. Wc first report results obtained us- 
ing the reaction field method for a spherical sample em- 
bedded in a long needle-shaped domain. We also report 
results fi'om simulations using the ES method for both 
a long needle-shaped sample and a spherical sample em- 
bedded in a long needle-shaped domain. 



low enough fields close to the classical phase transition, 
our perturbative Monte Carlo results, using the same re- 
action field method as in Ref. [l^, closely match the 
quantum Monte Carlo results from Ref. [l^. Using the 
reaction field method for = we get a Tc = 2.03 K, 
where Tc{B^ = 0) = fc{B^ = 0) since e{B^ = 0) = 1. 



1. Results from reaction field method 

To establish a comparison of the effective perturba- 
tive Hamiltonian with previous QMC resultsi^, we first 
performed Monte-Carlo simulations for a finite size sam- 
ple with open spherical boundary condition, containing 
N — 295 spins and with Jex in Eq. (Risindip) set to zero. 
These conditions are identical to the ones of Ref. p^ . 
As shown in Fig. [71 similarly to Ref. [11] , we used the 
reaction field criterion, set by the divergence of x when 
Xsph = ^ (see Eq. ((28)) ). to find the effective critical 
temperature Tc{Bx) as a function of the effective field 
Bx, where T and B^ are defined in Eqs. (fT5)l and (fTB)) . 
Xsph is calculated using 

Xsph = -^^ (Ml) , (33) 

where the prefactor a is given by 

a = 4-(5LMBC,,(0))' . (34) 

In the perturbative MC method, for determining (M|), 
we used the pseudo-operator defined by Eq. (|C4p . 
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FIG. 7: Finding Tc using the perturbative Monte-Carlo for 
a sphere of N=295 spins and Jcx = 0, by using the reaction 
field Xsph ~ ^ criterion at criticality 

The phase diagram as a function of the effective tem- 
perature T and the effective field B^, using the effective 
perturbative Hamiltonian (j22p and the above cumulant 
expansion is shown in Fig. [51 It can be seen that at 



2. Results from Ewald summation method — needle-shaped 
sample 

The simulations using the Ewald summation (ES) 
method were performed with simulation boxes of size 
L — 7,8, 9, with each box containing iV = 4 x L"^ spins. 
The dipolar interactions of ions inside the simulation 
boxes were derived via the ES technique and assuming 
an infinitely long needle-shaped sample^. That is, the 
additional demagnetization term correction from Eq. (j32l 
is not incorporated into the Ewald representation of the 
dipolar interactions between ions i and j. We determined 
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r(K) 

FIG. 8: Comparing the phase diagram of the perturbative 
Monte-Carlo with Quantum Monte Carlo results^^ as a func- 
tion of effective temperature and effective magnetic field for 
a sphere of N=295 spins and Jcx = , using the reaction field 
method of Ref. 

the critical temperature by finding the temperature at 
which the magnetization Binder ratio^, 

Q = l-l{Mt)/{M!)' , 

for system sizes L = 7,8, and 9 intersect. The intersec- 
tion point shown in Fig. [9l is at Tc ~ 1.92 K which is 
the critical temperature. (A/^) and (Af|) are calculated 
using Eqs. (|C4p and (|C9p within the perturbative effec- 
tive Hamiltonian scheme. As demonstrated in the inset 
of Fig. [H plotting Q as a function of L^/^iT - Tc) shows 
a good data collapse for system sizes L ~ 7,8, and 9, 
with the mean field exponent = 1/2. This is consistent 
with the argument that the upper critical dimension for 
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dipolar interactions is d = 3. A more rigorous analysis 
of three dimensional dipolar systems shows logarithmic 
finite size scaling corrections^^^. We have not investi- 
gated these corrections in this study as it is outside the 
scope of this work. As long as Tc{Bx 7^ 0) > 0, the crit- 
ical behavior should be controlled by the same classical 
critical exponents as for = 0. 



3. Results from Ewald summation method — spherical 
sample 

We have repeated the pcrturbative MC simulations us- 
ing the ES technique but with a slightly different twist 
to it. Instead of simulating a long needle-shaped bulk 
and using the Binder method to obtain the critical tem- 
perature, we simulate a sample with a spherical domain. 
We derived the effective dipolar interactions between the 
spins by using the ES technique for a spherical cavity 
The effect of the spherical boundary is taken into ac- 
count by incorporating the additional effective interac- 
tion of Eq. ([32|^ between spins i and j. Now, one can 
assume that this sphere is embedded in a long-needle- 
shaped bulk. Therefore, by recalling the derivation of 
Eq. ^ from Eq. ([221, where an effective field Bf"^ is 
applied to the magnetic moments of the sphere, one can 
determine the macroscopic x of the bulk, by calculating 
Xsph via ES method for a spherical sample. The pro- 
cedure that we use here is similar to the procedure one 
above that employed the reaction field for a finite size 
system and which led to the phase diagram in Fig. [S] 
The difference between the ES technique within a spheri- 
cal boundary and the reaction field method implemented 
in Ref. [l9j is that instead of using an open spherical 
boundary condition, and considering only bare dipolar 
interaction between a finite number of spins within a 
cutoff sphere, a simulation box with periodic boundary 
condition is considered. The effective dipolar interac- 
tions of ions inside the simulation box is derived via the 
Ewald summation technique. In this approach a spheri- 
cal boundary is considered for the whole simulation box 
and all the replicated images of the real box. In this case, 
each effective pairwise dipolar interactions described by 
the ES representation has added to it the extra inter- 
action term given by Eq. [32l Once again, the origin of 
this additional interaction is the demagnetization field, 
due to the polarization of the magnetic moments on the 
spherical boundary. In this approach, the system behaves 
much more like a macroscopic sphere compared to the one 
above that used the reaction field method. It is further 
assumed that this macroscopic sphere is embedded inside 
a macroscopic long macroscopic needle-shaped domain. 
Therefore, by employing the pcrturbative Monte Carlo 
method and using Eq. ([M)) . we Ccilculcite Xsph 

to obtain 

the critical temperature. Based on Eq. ((28)) . the critical 
temperature is calculated by finding where the Xsph = ^ 
criticality criterion is satisfied. As shown in Fig.O), for a 
simulation box oi L — 7 , we obtain Tr ~ 1.92 K for a zero 
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FIG. 9: (a) The Binder ratio crossing for L — 7, 8, 9 sys- 
tem sizes, performing MC and using ES technique for a long 
needle-shaped sample, with Bx ~ 0, Jgx = 0. T — T for 
Bx = 0- The inset shows that the Binder ratios collapse for 
the mean field exponent v — 1/2 to a very good degree, (b) 
Xsph calculated by performing MC simulation, using Eq. (|33p . 
The diamonds are for a finite size sphere using the reaction 
field scheme similar as in Ref. [l^ (i-^- same results as shown 
in Fig. [7] for the Bx = data). For the circles, we have ob- 
tained the interaction between the ions by the ES technique 
for L = 7 system size and incorporating the spherical bound- 
ary effect via the demagnetization term of Eq. (|32|l and using 
Bx = and Jox ~ 0, with again T — T for Bx = 0. As one 
can see the Tc ~ 1.92 K obtained here agrees with the Tc 
obtained using the Binder ratio crossing. 



transverse field and Jox = 0, very close to the Tc previ- 
ously derived using ES technique for a long needle-shaped 
sample and shown in Fig. [9^. Thus, the two approaches 
using ES technique lead to similar results. We believe 
that the difference between the classical Tc obtained via 
ES technique and the Tc{Bx = 0) obtained using the re- 
action field method — is because, in the reaction field 
method, the number of spins inside the cut-off sphere, 
which is embedded in the needle-shaped domain, is of 
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too limited size. By implementing Eq. ((28)) in the reac- 
tion field method, the effect of the spins on the spherical 
boundary for a limited size is in essence incorporated in 
a mean field manner in the simulation. For a limited 
size boundary, thermal fluctuations on the boundary are 
underestimated, hence resulting in an overestimated Tc- 
This overestimation of Tc, which decreases by increasing 
the size of the spherical boundary, is expected to vanish 
in the thermodynamic limit L s- oo. 



C. Nearest-Neighbor Exchange Interactions 

The zero transverse field critical temperature of 1.92 K 
obtained above lies quite far above the experimental criti- 
cal temperature of 1.53 K. As suggested by Chakraborty 
et at, it is reasonable to assume that the discrepancy 
may be related to a nearest-neighbor Heisenberg antifer- 
romagnetic exchange interaction. Indeed, in the related 
LiTbF4 material, it has long been known that a Jox cou- 
pling exists^. There has been no direct determination 
for the magnitude of this nearest-neighbor exchange in 
LiHoF4. However, there have been indirect estimations, 
considering Jox as a free parameter, such that the specific 
heat^ and susceptibility^ calculations based on mean 
field theory fit to the equivalent experimental measure- 
ments. Another procedure to determine Jcx, would be to 
fit theoretical calculation with neutron scattering data, 
similar to the procedure followed for LiTbF^^. Recently, 
R0nnow et. al^ have performed inelastic neutron scat- 
tering measurements on LiHoF4. Considering Jcx as a 
free parameter, they used the so called effective-medium 
theory to modify the mean field random phase approxi- 
mation parameters. They estimated Jox such that a best 
fit with the experimental phase diagram is obtained. For 
example, although for Jex=1.16 mK there is good agree- 
ment with experiment when 2.0 < < 4.0 Tesla, as 
is common in mean field theory calculations, the critical 
temperature is overestimated (by 14 percent) compared 
with the experimental critical temperature at zero ap- 
plied field B^ = 0. 

In our work here, we use Monte Carlo techniques and 
consider the exchange interaction as a free parameter. 
We can estimate the Jox strength by adjusting its value 
such that the experimental Tc is reproduced, as was 
done in Ref. [l^. Using the reaction field method per- 
formed for finite spheres in Ref. [l^, for N = 295 spins, 
Jox = 6.07 mK was obtained. As a check, we repeated 
our Monte Carlo simulations, also using the reaction field 
method for the same number of spins, and fitted Jox 
such that the experimental zero-field critical tempera- 
ture Tc = 1.53 K is reproduced. We obtained the same 
Jex = 6.07 mK as in Ref. [l^. It should be noted that, 
as reported in Ref. [l^ , one does not obtain a unique Jex 
value when performing simulations for different sphere 
sizes. The Jex value strongly depends on the number of 
spins considered. In Ref. (l9[, for the largest system size 
considered (N=3491), a Jox = 5.25 mK was required to 




r(K) 

FIG. 10: The Binder ratio crossing for L = 7, 8, 9 system 
sizes, performing MC and using ES technique for a cylindrical 
boundary with Bx ~ 0. Jox ~ 6.07 mK is set such that the 
critical temperature Tc « 1.53 K is obtained. T = T for 
Bx = 0. In the inset Xspii is calculated by performing Monte 
Carlo simulations, using Eq. (|33|l . The interaction between 
the ions is obtained by the ES technique for L — 7 system 
size and using a spherical boundary condition for B^ ~ 0. 
The same Jcx = 3.91 mK used and a similar Tc ~ 1.53 K is 
obtained. T = T for Bx = 0. 

obtain a Monte Carlo estimate of Tc of 1.53 K. There 
are two sources of errors that are affecting the value of 
the estimated Jox obtained by the reaction field method 
of Ref. [l^ . Firstly, for a given number of spins, when 
Monte Carlo simulations are performed to calculate Tc, 
the reaction field method estimates a higher value for 
Tc compared to the ES method. The sources of these 
errors are finite size effects and the underestimation of 
thermal fluctuations at the boundary, as we now explain. 
To push down the value of Tc obtained for Jox = such 
that it matches the experimental value for Tc, an anti- 
ferromagnetic Jox is required. For Jox = 0, the reaction 
field method generates a higher Tc compared to the ES 
method. Therefore, in order to push down the Tc ob- 
tained from Monte Carlo simulation to match the exper- 
imental value for Tc, a larger value for the antiferromag- 
netic Jex is required than the one required when using 
the ES method. Secondly, there is another source of er- 
ror affecting the value of the estimated Jex obtained by 
the reaction field method. It comes from the number of 
surface bonds, which depends on the radius of the cho- 
sen cut-off sphere. For ions close to the surface, some of 
the nearest-neighbors fall inside the spherical boundary 
while some remain outside. Because of the missing num- 
ber of exchange interactions on the boundary, the overall 
exchange estimated is forced to be larger than the actual 
value. When the ES technique is used in conjunction with 
periodic boundary conditions, this boundary effect prob- 
lem no longer exists, making the ES technique a more 
reliable tool for estimating Jc x^^''*^ - To estimate Jex us- 
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FIG. 11: The Binder cumulant crossing for L = 7, 8, 9 system 
sizes, performing perturbative MC and using ES technique 
for a long needle-shaped sample with Jcx = 3.91 mK. In (a) 
we have = 0.05 T and in (b) we have = 0.15 T. 



D. Transverse Field vs Temperature Phase 
Diagram 

Having determined a seemingly consistent value for 
Jex, we are now ready to perform Monte Carfo simu- 
lation for small transverse magnetic fields ■ The effect 
of quantum perturbations are incorporated through the 
effective Hamiltonian Eq. ([22|) . which is derived from the 
Sa;— rcscaled Hamiltonian Eq. p^ . To obtain the real 
temperature T and external transverse magnetic field 
from the effective values T and B^ used in the simulations 
we employ relations Eqs. and (|16|) . To illustrate the 
procedure, wc show the crossing of the Binder ratio Q for 
B.;^ = 0.05 T and B^ 0.15 T in Fig. [IH 

Interestingly, using each of the numerical methods dis- 
cussed above to obtain the phase diagram, it seems that 
for small B^ the final phase diagrams demonstrating the 
critical transverse field as a function of temperature are 
affected very little in respect to which specific technique 
is used. Figure [T^] shows the phase diagrams, using 
the perturbative Monte Carlo method implementing the 
reaction-field method and the Ewald summation tech- 
nique, compared with QMGi^ results and experiment^. 
We use Eq. ([T5]) and Eq. to obtain the real phys- 
ical transverse magnetic field, B^ and temperature T 
from T and B^- As one can see, all the phase diagrams 
obtained from the effective perturbative method show a 
good agreement with the quantum Monte Carlo result of 
Chakraborty et ali^, for small transverse fields up to a 
"real" physical transverse magnetic field B^ ~ 1.5 Tcsla, 
where we presume the lowest order cumulant formulation 
of the effective classical Hamiltonian model breaks down. 
This is the main result of this work. 

In conclusion, we confirm the results of Ref. [l^ but, 
perhaps unfortunately, we fail to explain the discrepancy 
between numerical and experimental results. We are thus 
led to ponder on theoretical reasons that may explain this 
discrepancy. We explore one such possibility in the next 
subsection and which is also the one that was put forward 
in Ref. [l^. 



ing our Monte Carlo simulations, we used the Binder ra- 
tio crossing method and employed both the ES technique 
for a long needle-shaped sample and the ES technique for 
a macroscopic sphere embedded in a long needle-shaped 
sample. For the latter case, the interactions of Eq. ([5^ . 
originating from the magnetic polarizations of the mag- 
netic moments on the spherical boundary were considered 
as well. The two Jex values so determined are the same, 
which is approximately Jcx = 3.91 mK, as illustrated in 
Fig. Uni Note that this value of Jcx = 3.91 mK is consis- 
tent with the one recently determined in Ref. [3J|. The 
definition of the exchange constant of 0.12 K in Ref. [s^l 
for Ising spins corresponds to JexCzz{Bx ~ 0) in our 
case. Using J^x ~ 3.91 mK and Czz{Bx = 0) = 5.51 



from Fig. [H we have JcxCzz (Bx 
ccllent agreement with Ref. [3J] 



= or 



0.119 K, in ex- 



E. Other Crystal Field Parameters 

As reported in Ref. [l^, we find that the numerical 
phase diagrams show a discrepancy with the experimen- 
tal phase diagram, even at asymptotically small trans- 
verse fields. Indeed, this was one of the main motiva- 
tions for the present work. As can be seen in Fig. 1121 our 
efforts in considering (i) a different Monte Carlo scheme 
and (ii) other ways to handle the long-range dipole-dipole 
interactions have not allowed us to resolve the discrep- 
ancy between the results from numerical simulations of 
Ref. [l^ and the experimental phase diagram of Ref. [l^ . 
Chakraborty et alr^ suggested that this discrepancy may 
be related to uncertainties in the crystal field parameters. 
We now briefly explore this possibility. 
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FIG. 12; The phase diagram of the critical transverse field as a function of temperature for LiHoF4. The closed boxes are the 
experimental phase digrarr^. The closed triangles are the phase diagram obtained by QMC— using the RF method for a finite 
sphere with A'^ — 295 spins. The open stars are the result from perturbative Monte Carlo (PMC) using the same reaction field 
(RF) method used in Ref. [l^ for a sphere with A'' = 295 spins. Quite importantly, as discussed in the text, the reaction field 
method leads to a considerable overestimate of Jox. The open circles are obtained, using the perturbative Monte Carlo in a 
needle-shaped domain using ES method. The open diamonds are obtained, using perturbative Monte Carlo in a bulk sphere 
embedded in a needle-shaped domain, using ES method and the spherical boundary effect of Eq. (|32p . 



As discussed in Appendix \^ crystal field parame- 
ters are usually obtained such that theoretical calcula- 
tions match with experimental data from electron para- 
magnetic resonance (EPR)^, inelastic neutron scattering 
(INS)^ or susceptibility measurements^^. Recalling the 
discussion that led to the derivation of the effective spin- 
1/2 description of LiHoF4 in Eq. ([TT|) . one realizes that 
the parameters Czz{Bx) and lS.{Bx) are implicit functions 
of the crystal field level energies and crystal field level 
wave functions. As a result, the mapping of the prob- 
lem to a spin-1/2 model depends on the chosen values of 
the i?" crystal field parameters (See Appendix |^ enter- 
ing in the description of the crystal field Hamiltonian Vc- 
This state of affairs is rendered particularly important, 
since, unfortunately, there appears to be some ambigu- 
ity in the literature about the empirical values of the B'^ 
parameters. All the numerical results that we obtained 
in the previous sections are based on the set of recent 
crystal field parameters obtained reported in Refs. p^ . 
which were also used in the stochastic series expansion 
quantum Monte Carlo of Ref. [l§| , and which were de- 
termined by fitting theoretically determined crystal field 
levels with those resolved in inelastic neutron scattering 
data. Recently, new electron paramagnetic resonance 
(EPR) spectroscopy experiments have been performed, 
in which the crystal field parameters were determined^. 
Based on the EPR data reported in Ref. [13], spectral 
parameters were refined in order to fit the observed de- 



pendencies of the resonance frequencies on the external 
magnetic field, giving a new set of crystal field parame- 
ters and an effective Lande g-factor reduced from its 
pure ^/g 5l ~ 5/4 value down to an effective — 1.21. 
Using this new set of crystal-field parameters, we ob- 
tain a different renormalization factor e{Bx) (Eq. ([T^ ) 
and effective transverse field Bx (Eq. P^ ) and as a re- 
sult, different Czz{Bx) and A(i3.^). One of the conse- 
quences of obtaining a different Czz , with the new CFP 
is that referring, to Eq. pT|) . a different = 0, Tc is 
obtained. Having determined a different Tc via this new 
set of CFP, the value of Jex required to match the ex- 
perimental Tc = 1.53 K is therefore different. In order to 
scrutinize "only" the effect of using a new set of CFP and 
to compare the phase diagram obtained using this new set 
of parameters with the results of Ref. [l9j in a rather sim- 
ple way, we repeated the perturbative Monte Carlo simu- 
lations, using the same reaction field method used above 
and done in Ref. [l^ for a finite size sphere of = 295 
spins and a newly determined Jox = 4.38 mK. At the end, 
after in essence repeating all the work discussed in Sec- 
tion [IVBTI a new phase diagram is derived. This phase 
diagram is plotted in Fig. [TH As it can be seen, this new 
phase diagram is consistent with the previous theoreti- 
cal work, {e.g. Ref. [l^ and Fig. [T2)) . Interestingly it 
therefore does not appear at this time that the different 
crystal field Hamiltonians available for LiHoF 4^^i^^i^^i'^^ 
are able to explain the significant discrepancy between 
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FIG. 13; Comparing Czz and A as a function of B^c, calcu- 
lated using two different crystal field parameters (CFP). The 
solid line is obtained, using Refs. [lill| CFP based on inelas- 
tic neutron scattering (INS) experiment. The dashed lines is 
obtained, using Ref. [33] based on electron paramagnetic res- 
onance (EPR) experiment. 
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FIG. 14: Comparing the phase diagrams of the critical trans- 
verse field as a function of temperature for LiHoF4, based 
on two different set of crystal field parameters. The closed 
triangles are the QMC results of Ref. [l^, using the RF 
method for a finite sphere with A'^ = 295, based on the CFP 
of Ref. [l^ . The open diamonds are obtained from our per- 
turbative Monte-Carlo, using the same RF method used in 
Ref. [3 for ^ sphere with A'^ — 295 spins, based on the CFP 
reported in Ref. [37| . 



the Bx — T phase diagram obtained by simulations com- 
pared to experimental results of Ref. [l3|. Finally, it 
should be emphasized that there is no difference in the 
results for this new set of CFP provided Jex is adjusted 
as well. On the other hand, different CFP lead to a sys- 
tematically different Tc if Jex is not adjusted. 



V. CONCLUSION 

With a perturbative Hamiltonian derived from a low 
energy effective spin- 1/2 description of LiHoF4, we have 
re-investigated the — T phase diagram with an inde- 
pendent approach for small B^/Tc where quantum fluc- 
tuations are weak. The method we used to incorporate 
pcrturbatively weak quantum fluctuations within a semi- 
classical Hamiltonian, because of its simple numerically 
tractable form, allows one to directly address possible fac- 
tors behind the discrepancy between results from experi- 
ments and from classical Monte Carlo simulations in the 
vicinity of Tc- This method can be easily generalized to 
more complicated quantum magnetic Ising models, where 
the Ising-like term is the dominant term and the other 
non commuting terms are considered as weak perturba- 
tions. In particular, if one is interested in studying nu- 
merically the effect of nonzero Bx in the diluted regime 
of LiHoa;Yi_2,F4, this perturbative method should be di- 
rectly applicable by performing Monte Carlo simulations 
of the appropriate low energy Hamiltonia n^^i^^ . 

To perform semi-classical Monte Carlo simulations 
that handle the magnetostatic long-range dipole-dipole 
interactions properly, we applied the Ewald summation 
technique for two different geometries. In order to deter- 
mine Tc, we used the Binder magnetization ratio cross- 
ing for a long needle-shape sample, and we used the 
Xsph = ^ criterion for a spherical sample embedded in- 
side a long needle-shaped domain. We obtained the same 
Tc for both cases and, consequently, determined the same 
value for Jex. The values of the Tc and Jex that we calcu- 
lated are somewhat different from the Tc and Jex values 
found in Ref. [l^. This difference originates from using 
open boundary conditions and a finite spherical cutoff 
in Ref. [T^, which underestimates the thermal fluctua- 
tions at the boundary. We found that although we used 
a different method and found a different Jqx, the flnal 
Bx — T phase diagram obtained here is the same in the 
low Bx/Tc limit as in the previous resultsi^. As a result, 
we tentatively conclude that the discrepancy between the 
theoretical and experimental results is not of computa- 
tional origin. To explore a possible explanation for the 
discrepancy, we considered a different set of crystal field 
parameters. 

A consideration of different crystal field parameters 
(CFP), which lead to a different estimate for Jex does 
not, however, at the end produce a dramatically differ- 
ent Tc vs Bx phase diagram. This preliminary result that 
only considers one set of alternative CFP goes against 
the suggestion of Ref. [l^ that a possible origin of the 
discrepancy might be due to the ambiguity in CFP. It 
is perhaps surprising that the consideration of a rather 
different set of CFP compared to those used in Ref. [Toj 
affects the phase diagram so little once Jex has been re- 
adjusted to match the experimental Tc{Bx = 0) = 1.53 K 
value. Therefore the origin of the discrepancy between 
numerics and experiment remains fully unexplained. 

The method we obtained in the present work could 
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be used to carry on further investigation of the cause 
of the discrepancy. Without this tool, it would have 
been somewhat less straightforward to have investigated 
the relevance of the various factors that we investigated 
in this paper. The disagreement with the experimental 
phase diagram of Rcf. would suggest that it may 

be worthwhile to revisit the experimental determination 
of the Bx vs Tc phase diagram. On the other hand, in 
both the work presented here and in that of Ref. p^ . 
a very simple spin Hamiltonian was considered. Specif- 
ically, only long-range magnetostatic dipole-dipole and 
isotropic (Heisenberg) nearest-neighbor exchange inter- 
actions were considered. The faster decreasing Tc{Bx), 
compared to the experimental case indicates that per- 
haps there are effects at play in the real material that 
weaken quantum fluctuations for small B^- 

In other words, there may be other couplings in the ef- 
fective theory in addition to those in the simplest trans- 
verse field Ising model (TFIM) of Eq. ([TT]). As illustrated 
in Fig. O the terms that we ignored when passing from 
Eq. (jTO]) to Eq. pT|) seem too small to be able to re- 
solve this issue. It might be necessary to consider the 
possibility that not completely negligible anisotropic ex- 
change, higher order multipolar exchange interaction, or 
magneto-elastic couplings may be at play in LiHoF4. 

Finally, wc note that it would be interesting if one 
could study other magnetic materials similar to the 
LiIIoF4 compound and that could provide another re- 
alization of a TFIM. Recently, a mean-field theory cal- 
culation has concluded that Ho(OII)3, which is an insu- 
lating hexagonal dipolar Ising ferromagnet, is very well 
described by a TFIM when a magnetic field B^ is applied 
perpendicular to the Ising spin direction^. This mate- 
rial constitutes a close analogue of LiHoF4 and, when 
diamagnetically diluted with , may potentially be an 
analogue of LiHoa;Yi_a;F4. The existence of another ex- 
perimental candidate for the study of the TFIM with 
long-range dipolar interaction presents the opportunity 
to re-investigate the puzzling properties of pure and di- 
luted LiIIoF4 in a new material, shedding light on the 
physics of dipolar Ising systems in both zero and nonzero 
applied transverse field. The method we have employed 
in this work is a suitable tool to study these new pro- 
posed quantum magnetic Ising materials beyond mean 
field theory and provides a tool to make comparison with 
future experiments performed on these proposed TFIM 
materials. 

To conclude, we hope that the work presented here 
stimulates further theoretical and experimental studies of 
LiIIoF4 in the regime of small transverse field B^ where 
the classical paramagnetic to ferromagnetic transition is 
only perturbatively affected by B^- 
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APPENDIX A 

In this Appendix we briefly discuss how the crystal 
field Hamiltonian of LiIIoF4 is written in terms of angular 
momentum operators and crystal field parameters. 

In the point charge approximation description of the 
crystal field, we assume that the ions interacting with 
Ho'^"'" electrostatically are close to point charges. The 
potential at r is simply the sum of point charge coulomb 
interaction potential 



(Al) 



where is the position and the total electric charge of 
the z'th ion. V{t) can be expanded as 



nr,0,</>) = ^^r"7„ 



n— a 



where 



E 



no. \yi 1 H^i ) 



(2n+l) i?; 



n+l 



(A2) 



(A3) 



and the Z„q,'s are the spherical harmonics containing 
sin0 or cos</<^. To get the crystal field Hamiltonian Vc^ 
one must sum this energy over all of the valence electrons 
of the holmium (Ho'^+) moments, hence we have: 



Vr. 



(A4) 



According to arguments provided by Stevens^S for evalu- 
ating the matrix elements of the crystal field Hamiltonian 
between wave functions specified by the angular momen- 
tum J, the crystal field Hamiltonian can be written in 
term of Stevens' operator equivalents O", built out of 
the vector components of J operators. 



Vc = Y^BlOl . 



(A5) 



The Stevens' equivalent operators act on the angular mo- 
mentum states of the wave functions. The matrix ele- 
ment of the radial part of the wave function is incorpo- 
rated in the parameters, usually determined by fitting 
to experimental (e.g. spectroscopic) dat a^^i^^'^^ . From 
angular momentum algebra, in the case of 4/ electrons, 
we need to consider only n = 0, 2, 4, 6 in the sum (jASp . 
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The choice of coefficients in Hamihonian (|A5p that 
do not vanish and have nonzero corresponding matrix el- 
ements is dictated by the point symmetry group of the 
crystalline environment. The details of the method and 
conventions for expressing the crystal field Hamiltonian 
can be found in the review paper by Hutchinga^. The 
point group symmetry of LiHoF4 has ^4 symmetry, which 
means the lattice is invariant respect to a ^ rotation 
about the z axis and reflection with respect to the x — y 
plane. The crystal field Hamiltonian for LiHoF4 is there- 
fore written as 

Vc = BlOl+BlOl + BfOf +BfOf 

+BlOl+Bfor + Bt'or. (A6) 

The relevant _B" crystal field parameters must be deter- 
mined experimentally The relevant operator equivalents 
are given in terms of angular momentum operators^ (J^, 
J+, J_, J2) by 

0° = ij'l - j2 , 

O4 = 35 Jz"* - 30 jVf -f- 2bJl - QJ^ -f iJ^ , 

Of = \{A + Jt), 
Of = \{A-Jt). 

Ol = 231J| - 315 JV^^ + 735 + 105 JV^2 

-525 J^J^ + 294J2 _ 5 j6 + 40 j4 _ 60 , 

Of = + Jl)(llj2- J2_38) + H.c. ,and 

Of = ^(J|- Jl)(llJ,2- j2-38) + H.c. (A7) 

The parameters are chosen such that the resulting 
energy levels match those determined from spectroscopic 
data. Two different set of experimentally determined 
crystal field parameters are given in Table HI The first 
set of the parameters was determined by inelastic neu- 
tron scattering reported in Ref. and implemented in 
the calculations of Ref. [l^. The next set of B" pa- 
rameters were determined using electron paramagnetic 



resonance (EPR) spectroscopy, and reported in a recent 
work^. 



APPENDIX B 

In this Appendix, starting from Eq. ()2ip . we give 
the details of the derivation of the effective perturba- 
tive Hamiltonian Heg [tpi] by cumulant expansion, when 
quantum fluctuations are small. Deriving H^s [ipi] , as de- 
fined by Eq. (|2Gp one can rewrite the partition function 
of the system in a classical form. 

Referring to Eq. (|21|). recalling that \^), is a direct 
product of erf eigenstates, the expectation value {iplo'x IV') 
is zero, so (iplHlip) = (-iAI^^olV')- Defining Eo{ip) = 
(V'l-ffolV'), we can write (?M (W - (V'lWlV))" |V'> = 



Parameter 


Ref. [15] Ref. [37] 




-0.696 K -0.609 K 


Bl 


4.06 X 10"^ K 3.75 x 10"^ K 


Bf 


4.18 X 10"^ K 3.15 X 10"^ K 


Bf 


K 2.72 X 10"^ K 


Bi 


4.64 X 10"® K 6.05 x 10"® K 


Bf 


8.12 X 10"'' K 6.78 x 10"'' K 


Bf 


1.137 X 10"" K 4.14 X 10"" K 



TABLE I: The first column is the crystal field parameters 
(CFP) for LiHoF4 determined experimentally by fitting the 
results of random phase approximation spin-wave dynamics 
calculation to neutron scattering data from Ref. [15|. The 
second column is the crystal field parameters estimated using 
electron paramagnetic resonance (EPR) spectroscopy experi- 
menli^. 



Performing a polynomial expansion on (7i — i?o('0))" = 
[{Hq — Eq{i(,')) + Hi]"^ , and keeping terms to order of 
O(F^) in the polynomial expansion {Hi oc F), we have 



rii .712 ■'>^3 

= {iP\Hi [Ho - So(V)]""' i^ilV^) . (Bl) 

I 

To write Eq. (jBip we have used the fact that and 

mHo - EoWrHiiHo - Eom^li^) = , (B3) 

for integer numbers m and k. The effect of af on \tjj) is 
(?M (Ho - Eoi^)y' IV') = (B2) to flip the spin i. We define erf IV") = where f.i/j 
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means that the i'th spin has flipped, such that if the i'th 
spin was in the 1 1) or the 1 1) eigenstate of erf, it changes 
into the \ l) or | f) state respectively. In Eq. (jBip . using 
Hi = -rX; .CTf , we get 

= E^,,{f^i'\ [Ho EoWr' \m ■ (B4) 

Here {fii/j\ [Hq — i?o(^/')]"~^ l/j^) is zero, unless i ~ j. 



Thus, Eq. (|B1[) can be written as 

{H - i?o(V))" IV') = ^ [EoiM) - EoWr-^ -(BS) 

i 

Considering the definition of HeS, by substituting 
EoiM) - Eoiij) = 2{hi + ho)af in Eq. 1^, we obtain 

OO ^ 

= ifo-/3r2^^-[-2/3(/^, + /io)]""' 

i n> 1 

= Fo + ^1 [2/3(/i« + ho)] 

i 

~Fo [2(3{h, + ho)]}. (B6) 

In Eq. (|B6[) . /i^ is the total local field affecting the spin 
at site i by other spins , which is 

/ij - ^ C^jCTj - Jox ^ ctnn , (B7) 

j^i NN 

and ho is the external longitudinal field in the z direction. 
The functions Fo{x) and Fi{x) are defined as 



Foix) = 
Fi{x) 



cosh(a;) — 1 

X'^ 

sinh(a;) — x 



(B8) 



APPENDIX C 

In this Appendix, we establish the relationship be- 
tween the real thermodynamical quantities as physical 
observables and their corresponding pseudo-operators, 
which are obtained using the pcrturbative effective clas- 
sical Hamiltonian of Eq. ((22)) . These thermodynami- 
cal quantities are calculated by employing the derived 
pseudo-operators in our perturbative classical Monte 
Carlo simulations. 

Writing the partition function in terms of the perturba- 
tive effective Hamiltonian iJoff [V't]; pseudo-operators 
corresponding to (E), (A-/^), (A4), (M^), and (A/^), 
which should be calculated to obtain thermodynamical 



quantities using Monte Carlo simulations are written as 

(CI) 

(C2) 

- (-^) (C3) 



{Ml 



dV 



dho J P dhl 



4/3^ 



dho 



(C4) 
(C5) 



-6/3' 



dh?. \ dho 



dh?, 



off 



4 / dHcs 



dho 



(C6) 



where E, Mz and Mx are the energy and magnetization 
in the z and x direction and their equivalent pseudo- 
operators which should be calculated are on right. The 
thermal average is denoted by (...) . Applying the deriva- 
tives and using the perturbative effective Hamiltonian 
dm), we find: 

E = ^0 + 2/3r2 ^KFi [2/3(/i, + M] 

i 

^FoW{h, + ho)]} 

-f/jr^ 2/?(/i« + ho){cjtF[^^ [2/3(/i, + ho)] 

i 

-F^'^ + ho)]} , (C7) 

while Mx is 

Mx = ~2prJ2W!Fi[2pih,+ho)] 

i 

-Fo [2p{h, + ho)]} (C8) 

and Mz 

Mz =T.^l^ 2/32r2 Y.{atF['^ mK + ho)] 

i i 

[2f3{h, + ho)]} , (C9) 

with f}"\x) defined as f/"^ = where i = 1, 0. 

In order to find an expression for (M^) and {M^), we 
need to calculate ^5^, and . We 

find: 

dho 



^^o^^) [2/3(/j,; + ho)]} , 



-i^o*'^ [2/3(/j, + ho)]} , 



(CIO) 



(Cll) 
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dhl 



i^o<^^ [2/3(/i, + /lo)]} and (C12) 



16/35r2^KFf) [2/3(/i. + M] 
-i^o<") [2/3(/i, + ho)]} . (C13) 
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